% default_time_sim

function [tau,Tot_out_cap_loan,Tot_periodic_loss]=default_time_sim_auto(T,ZC,Recov,Lambda,ref_ent,k,Rho)

% clear;
% clc;


% T=30; % Maturity of the loan or repayment (load T=4, repayment T=30)
% ZC=5; % Constant Interest Rate
% Recov=40; % Recovery rate
% Lambda=1;  % Number of failure per year

% ref_ent=100; % Nb of obligors
% k=1000; % Nb of simulations
% Rho=0.3; % Correlation
Rho=Rho/100;
ZC=ZC/100; % Constant Interest Rate
% % ZC=0.05; % Constant Interest Rate
Recov=Recov/100; % Recovery rate
% % Recov=0.4; % Recovery rate
Lambda=Lambda/100; % Number of failure per ye
% % Lambda=.02; % Number of failure per year
Amount=zeros(ref_ent,1); % Vector of notional amount for each credit


Amount(:)=100; %each credit has a notional amount of 100 units



% Hazard Rate - 5 nods

xx=min((Lambda-.005)/2,(.995-Lambda)/2);
delta=min(xx,0.02); %Hazard Rate Cycle
hazard=(max(Lambda-2*delta,.005):delta:min(Lambda+2*delta,.995)); %make the hazard rate in the center and extend to 5 points

% correlation of the individual default
% delta2=min(min((Rho-.05)/2,(.95-Rho)/2),.05); %Correlation cycle between .05 and .95 with 5 nods
% R=(max(Rho-2*delta,.05):delta2:min(Rho+2*delta,.95));
R=Rho;
mu=zeros(ref_ent,1);
% for R_cycle=1:5 % For each correlation cycle
for R_cycle=1 % For each correlation cycle
    % Compute teh correlation matrix
    corr=ones(ref_ent,ref_ent)*R(R_cycle) + (1-R(R_cycle))*eye(ref_ent);
    %************************************************************************
    % Generate default times with Gaussian Copula and Intensity
    rndmat1=(mvnrnd(mu,corr,k));
    def_t=-log(normcdf(rndmat1));
    %*********************************************************************
    S_fees=zeros(5,3);%Dummy for memorizing the simulated payment leg
    S_default=zeros(5,3);%Dummy for memorizing the simulated default leg
    M_fees=zeros(5,3);%variable which memorize the payment leg for each
    % loop of recovery and correlation
    M_default=zeros(5,3);%variable which memorize the payment lef for each
    % loop of recovery and correlation
    Total_loss=zeros(k,5); % Vector that keeps total losses for each simulation
    
    S_fees_loan=zeros(5,1);%Dummy for memorizing the simulated payment leg for the loan fund
    S_default_loan=zeros(5,1);%Dummy for memorizing the simulated default leg for the loan fund
    M_fees_loan=zeros(5,1); % Idem
    M_default_loan=zeros(5,1); % Idem
    
    for n=1:k % start the simulation loop
        %         for def_cycle=1:5 %Start the recovery loop
        for def_cycle=3 %Start the recovery loop
            [time,index]=sort(def_t(n,:)); % Sort Copula Gaussian Time for each simulation
            tau=time ./ hazard(def_cycle); % Generate a vector of default time adjusted for hazard rate
            
            % Calculate the total loss in the k-th simulation
            tot_loss=0;
            loss=zeros(ref_ent,1);
            for i=1:ref_ent
                % If the simulated default time for the generic credit is < than the
                % CDO maturity, there is a loss
                if tau(1,i) < T % For each time default, checks if it happen before T
                    loss(i)=(1-Recov)*Amount(1); % Each Obligator has the same capital and the same recovery rate
                    tot_loss=tot_loss+loss(i); % Sum the individual losses
                end
            end
            Total_loss(n,def_cycle)=tot_loss;
            
            %************************************************************************
            % Loan Fund
            %************************************************************************
            % Discount leg
            
            PV_def=0;
            PV_premium=0;
            periodic_loss=zeros(T,1);%Stores the accumulated loss at each payment date
            out_capital=zeros(T,1);%Outstanding tranche capital at each payment date
            fee=zeros(T,1);
            total_fee=0;
            indicator=0;
            c=0;
            if tot_loss <= 0
                PV_def = 0;
            else
                for i=1:ref_ent
                    if tau(1,i) < T
                        indicator=indicator+loss(i); % We memorize the cumulative losses
                        if c==0; % First Loss
                            r=ZC; % Discount factor at default
                            disc_fact_def=(1+r)^(-tau(1,i));
                            PV_def=PV_def+indicator*disc_fact_def;
                            c=1;
                        else
                            r=ZC;
                            disc_fact_def=(1+r)^(-tau(1,i));
                            PV_def=PV_def+loss(i)*disc_fact_def;
                        end
                    end
                end
            end
            
            % Premium Leg Simulation
            
            for i=1:T
                periodic_loss(i)=0;
                for j=1:ref_ent
                    if tau(1,j) < i;
                        % Calculate the accumulated portfolio losses
                        periodic_loss(i)=periodic_loss(i)+(1-Recov)*Amount(1);
                    end
                end
                % Outstanding capital at each payment date
                out_capital(i)=max(sum(Amount)-periodic_loss(i),0);
                fee(i)=((1+ZC)^(-i))*out_capital(i);
                PV_premium=PV_premium+fee(i); %DV01
            end
            
            
            %*************************************************************************
            if (R_cycle==1) && (def_cycle==3) % We compute only total periodic losses for the mid scenarios
                Tot_out_cap_loan(n,:)=out_capital';
                Tot_periodic_loss(n,:)=periodic_loss';
            end
            %**************************************************************************
            
        end
    end
end

% total fund left after 30 years (1000)simulation result
% hist(Tot_out_cap_loan(1:1000,30));figure(gcf);
% figure;
% 
% hist(ref_ent*Amount(1)-Tot_out_cap_loan(1:1000,30))*100/(ref_ent*Amount(1));